#include <stdio.h>
#include <math.h>
#include <curand.h>
#include <curand_kernel.h>

#include "heap.cuh"
#include "bench.cuh"


/********************************************************
*	CUDA optimized SSA, minimized model with competition;
*	source ssatest.c
*********************************************************/
__global__ void ssaKernel6(curandState * STATE, heapElement * ELEMENT, int seed){

	//Register
	#include "exp6_parallel_imp\ssaspecies.h"

	//Random and indexing;
	int THINDEX = threadIdx.x;
	int BLINDEX = blockIdx.x;
	
	//Parameters
	float ta				= 1.0*ELEMENT[BLINDEX].vector[0];				//binding probability
	float tda				= 0.0001*ELEMENT[BLINDEX].vector[1];			//unbinding probability
	float tdg				= 0.00001*ELEMENT[BLINDEX].vector[2];			//degradation probability
	float r_vp16_vp16		= 0.0001*ELEMENT[BLINDEX].vector[3];			//fully activated transcription
	float r_vp16_krab		= 0.01*r_vp16_vp16;								//competitive binding
	float tk				= 4.0;											//KRAB right of way
	float rdg				= 0.00001*ELEMENT[BLINDEX].vector[4];			//reporter degradation
	float cta				= 0.5*r_vp16_vp16*0.1;							//constitutive promoter transcription, LOWER FREQUENCY

	//Float
	float fbfp = ELEMENT[BLINDEX].vector[5];
	float fmct = ELEMENT[BLINDEX].vector[6];
	
	//Sample points
	float t2	= 36000; float t4 = 72000; float t6 = 108000; float t8 = 144000; float t10 = STEPCOUNT; 
	float yb_s	= 0;
	float ym_s	= 0;
	
	//Induction
	int induced = 1 + 1*(THINDEX > HISTSIZE/2.0);

	//Target values
	float yb_t2 = (induced == 1)*0.003 + (induced == 2)*0.184;
	float yb_t4 = (induced == 1)*0.038 + (induced == 2)*0.101;
	float yb_t6 = (induced == 1)*0.050 + (induced == 2)*0.066;
	float yb_t8 = (induced == 1)*0.063 + (induced == 2)*0.023;
	float yb_t10 = (induced == 1)*0.110 + (induced == 2)*0.015;

	float ym_t2 = (induced == 1)*0.149 + (induced == 2)*0.005;
	float ym_t4 = (induced == 1)*0.093 + (induced == 2)*0.075;
	float ym_t6 = (induced == 1)*0.059 + (induced == 2)*0.078;
	float ym_t8 = (induced == 1)*0.031 + (induced == 2)*0.088;
	float ym_t10 = (induced == 1)*0.012 + (induced == 2)*0.273;
	


	//Load vectors
	#include "exp6_parallel_imp\ssacoefs.h"
	
	
	//High log
	#define HL_AGM 20
	#define HL_PI 3.14159
	#define HL_LN2 0.69314
	#define HL_MBITS 8
	float hl_f;
	float hl_a;
	float hl_g;
	float hl_at;
	float hl_s;
	float hl_x;
	float hl_y;
	float hl_lnf;

		
	//Misc
	//int i		= 0;
	int j		= 0;
	int mi		= 0;
	int tstep	= 0;
	float tau	= 0;
	float t	= 0;
	
	//Ssa misc
	float r1;
	float r2;
	float wsum = 0;
	#include "exp6_parallel_imp\ssawthreshinit.h"
	float wtcumsum = 0;

	//Random number generator;
	unsigned z1;
	unsigned z2;
	unsigned z3;
	unsigned z4;

	z1 = 128+seed*THINDEX*BLINDEX;
	z2 = 129+seed*THINDEX*BLINDEX;
	z3 = 130+seed*THINDEX*BLINDEX;
	z4 = seed*THINDEX*BLINDEX;

	while(tstep <= STEPCOUNT){

		r1 = HybridTaus(&z1,&z2,&z3,&z4);
		r2 = HybridTaus(&z1,&z2,&z3,&z4);

		//Include propensities calculations
		#include "exp6_parallel_imp\ssaweights.h"
	
		//Calculate wsum
		//weight thresholds
		wsum = 0;
		wtcumsum = 0;
		wsum = wsum  + w0 + w1 + w2 + w3 + w4 + w5 + w6 + w7 + w8 + w9 + w10 + w11 + w12 + w13 + w14 + w15 + w16 + w17 + w18 + w19 + w20 + w21 + w22 + w23 + w24 + w25 + w26 + w27 + w28 + w29 + w30 + w31 + w32 + w33 + w34 + w35 + w36 + w37 + w38 + w39 + w40 + w41;
		#include "exp6_parallel_imp\ssawthresh.h"


		////Determine tau;
		//tau = 1.0/wsum * logf(1.0/(1-r1));
		hl_f = 1/(1-r1);
		hl_s = hl_f*pow(2.0,HL_MBITS);		
		hl_y = 4.0/hl_s;
		hl_x = 1;
		//agm start
		hl_a = 0.5*(hl_x+hl_y);
		hl_g = sqrt(hl_x*hl_y);
		for(j=0; j<HL_AGM; j++){
			hl_at = hl_a;
			hl_a = 0.5*(hl_a+hl_g);
			hl_g = sqrt(hl_g*hl_at);
		}
		//agm end
		hl_lnf = HL_PI / (2*hl_a) - HL_MBITS * HL_LN2;
		tau = 1/wsum * hl_lnf;

		////Determine reaction index.
		mi = -1;
		//for (i = 0; i<REACTIONS; i++)	mi = (r2 < wthresh[i] && mi == -1)*i + (r2 >= wthresh[i] || mi != -1)*mi;
		#include "exp6_parallel_imp\ssawtmi.h"
		

		////Update species matrix;
		//for (i=0; i<SPECIES; i++) species[i] = species[i] + dif[mi][i];
		#include "exp6_parallel_imp\ssadifmx.h"
		

		//Update time and step
		t = t+tau;
		tstep++;

		//switching
		//1: er_pi
		//2: pi_er
		species26 = (induced == 1 && tstep > t2)*INDMAX + (induced == 2 && tstep < t2)*INDMAX;
		species27 = (induced == 2 && tstep > t2)*INDMAX + (induced == 1 && tstep < t2)*INDMAX;

		//Sample values
		yb_s += (tstep == t2)*((species0 * fbfp)-yb_t2)*((species0 * fbfp)-yb_t2);
		yb_s += (tstep == t4)*((species0 * fbfp)-yb_t4)*((species0 * fbfp)-yb_t4);
		yb_s += (tstep == t6)*((species0 * fbfp)-yb_t6)*((species0 * fbfp)-yb_t6);
		yb_s += (tstep == t8)*((species0 * fbfp)-yb_t8)*((species0 * fbfp)-yb_t8);
		yb_s += (tstep == t10)*((species0 * fbfp)-yb_t10)*((species0 * fbfp)-yb_t10);

		ym_s += (tstep == t2)*((species1 * fmct)-ym_t2)*((species1 * fmct)-ym_t2);
		ym_s += (tstep == t4)*((species1 * fmct)-ym_t4)*((species1 * fmct)-ym_t4);
		ym_s += (tstep == t6)*((species1 * fmct)-ym_t6)*((species1 * fmct)-ym_t6);
		ym_s += (tstep == t8)*((species1 * fmct)-ym_t8)*((species1 * fmct)-ym_t8);
		ym_s += (tstep == t10)*((species1 * fmct)-ym_t10)*((species1 * fmct)-ym_t10);


	}
	
	ELEMENT[BLINDEX].errorbfp[THINDEX] = yb_s;
	ELEMENT[BLINDEX].errormct[THINDEX] = ym_s;
	
	
}


__global__ void ssaKernel8(curandState * STATE, heapElement * ELEMENT, int seed){

	//Register
	#include "exp8_parallel_imp\ssaspecies.h"

	//Random and indexing;
	int THINDEX = threadIdx.x;
	int BLINDEX = blockIdx.x;
	
	//Parameters
	float ta				= 1.0*ELEMENT[BLINDEX].vector[0];					//binding probability
	float tda				= 0.0001*ELEMENT[BLINDEX].vector[1];				//unbinding probability
	float tb				= 1.0*ELEMENT[BLINDEX].vector[2];					//binding probability
	float tdb				= 0.0001*ELEMENT[BLINDEX].vector[3];				//unbinding probability
	float tdg				= 0.00001*ELEMENT[BLINDEX].vector[4];				//degradation probability
	float r_vp16_vp16		= 0.0001*ELEMENT[BLINDEX].vector[5];				//fully activated transcription
	float r_vp16_krab		= 0.01*r_vp16_vp16*ELEMENT[BLINDEX].vector[12];		//competitive binding
	//float tk				= 4.0*ELEMENT[BLINDEX].vector[11];					//KRAB right of way
	float tk				= 4.0*ELEMENT[BLINDEX].vector[14];												//KRAB right of way
	float rdg				= 0.00001*ELEMENT[BLINDEX].vector[6];				//reporter degradation
	float cta				= 0.5*r_vp16_vp16*0.1*ELEMENT[BLINDEX].vector[10];	//constitutive promoter transcription, LOWER FREQUENCY
	float p_leak_pir		= 1.0*0.1*ELEMENT[BLINDEX].vector[7];				//pir promotor leakage
	float p_leak_etr		= 1.0*0.1*ELEMENT[BLINDEX].vector[8];				//etr promotor leakage
	float tas				= 1.0*ELEMENT[BLINDEX].vector[11];					//TAL A/B asymmetry
	float pas				= 1.0*ELEMENT[BLINDEX].vector[13];					//Induction asymmetry
		

	//Float
	float fbfp = 0.001*ELEMENT[BLINDEX].vector[9];
	float fmct = 0.001*ELEMENT[BLINDEX].vector[9];
	
	//Sample points
	float t2	= 36000; float t4 = 72000; float t6 = 108000; float t8 = 144000; float t10 = STEPCOUNT; 
	float yb_s	= 0;
	float ym_s	= 0;
	
	//Induction
	int induced = 1 + 1*(THINDEX > HISTSIZE/2.0);

	//Target values
	float yb_t2 = (induced == 1)*0.003 + (induced == 2)*0.184;
	float yb_t4 = (induced == 1)*0.038 + (induced == 2)*0.101;
	float yb_t6 = (induced == 1)*0.050 + (induced == 2)*0.066;
	float yb_t8 = (induced == 1)*0.063 + (induced == 2)*0.023;
	float yb_t10 = (induced == 1)*0.110 + (induced == 2)*0.015;

	float ym_t2 = (induced == 1)*0.149 + (induced == 2)*0.005;
	float ym_t4 = (induced == 1)*0.093 + (induced == 2)*0.075;
	float ym_t6 = (induced == 1)*0.059 + (induced == 2)*0.078;
	float ym_t8 = (induced == 1)*0.031 + (induced == 2)*0.088;
	float ym_t10 = (induced == 1)*0.012 + (induced == 2)*0.273;
	


	//Load vectors
	#include "exp8_parallel_imp\ssacoefs.h"
	
	
	//High log
	#define HL_AGM 20
	#define HL_PI 3.14159
	#define HL_LN2 0.69314
	#define HL_MBITS 8
	float hl_f;
	float hl_a;
	float hl_g;
	float hl_at;
	float hl_s;
	float hl_x;
	float hl_y;
	float hl_lnf;

		
	//Misc
	//int i		= 0;
	int j		= 0;
	int mi		= 0;
	int tstep	= 0;
	float tau	= 0;
	float t	= 0;
	
	//Ssa misc
	float r1;
	float r2;
	float wsum = 0;
	#include "exp8_parallel_imp\ssawthreshinit.h"
	float wtcumsum = 0;

	//Random number generator;
	unsigned z1;
	unsigned z2;
	unsigned z3;
	unsigned z4;

	z1 = 128+seed*THINDEX*BLINDEX;
	z2 = 129+seed*THINDEX*BLINDEX;
	z3 = 130+seed*THINDEX*BLINDEX;
	z4 = seed*THINDEX*BLINDEX;

	while(tstep <= STEPCOUNT){

		r1 = HybridTaus(&z1,&z2,&z3,&z4);
		r2 = HybridTaus(&z1,&z2,&z3,&z4);

		//Include propensities calculations
		#include "exp8_parallel_imp\ssaweights.h"
	
		//Calculate wsum
		//weight thresholds
		wsum = 0;
		wtcumsum = 0;
		wsum = wsum  + w0 + w1 + w2 + w3 + w4 + w5 + w6 + w7 + w8 + w9 + w10 + w11 + w12 + w13 + w14 + w15 + w16 + w17 + w18 + w19 + w20 + w21 + w22 + w23 + w24 + w25 + w26 + w27 + w28 + w29 + w30 + w31 + w32 + w33 + w34 + w35 + w36 + w37 + w38 + w39 + w40 + w41;
		#include "exp8_parallel_imp\ssawthresh.h"


		////Determine tau;
		//tau = 1.0/wsum * logf(1.0/(1-r1));
		hl_f = 1/(1-r1);
		hl_s = hl_f*pow(2.0,HL_MBITS);		
		hl_y = 4.0/hl_s;
		hl_x = 1;
		//agm start
		hl_a = 0.5*(hl_x+hl_y);
		hl_g = sqrt(hl_x*hl_y);
		for(j=0; j<HL_AGM; j++){
			hl_at = hl_a;
			hl_a = 0.5*(hl_a+hl_g);
			hl_g = sqrt(hl_g*hl_at);
		}
		//agm end
		hl_lnf = HL_PI / (2*hl_a) - HL_MBITS * HL_LN2;
		tau = 1/wsum * hl_lnf;

		////Determine reaction index.
		mi = -1;
		#include "exp8_parallel_imp\ssawtmi.h"
		

		////Update species matrix;
		#include "exp8_parallel_imp\ssadifmx.h"
		

		//Update time and step
		t = t+tau;
		tstep++;

		//switching
		//1: er_pi
		//2: pi_er
		species26 = (induced == 1 && tstep > t2)*INDMAX + (induced == 2 && tstep < t2)*INDMAX;
		species27 = (induced == 2 && tstep > t2)*INDMAX + (induced == 1 && tstep < t2)*INDMAX;

		//Sample values
		yb_s += (float)pow((float)(tstep == t2)*((species0 * fbfp)-yb_t2)*((species0 * fbfp)-yb_t2), (float)0.5); 
		yb_s += (float)pow((float)(tstep == t4)*((species0 * fbfp)-yb_t4)*((species0 * fbfp)-yb_t4), (float)0.5);
		yb_s += (float)pow((float)(tstep == t6)*((species0 * fbfp)-yb_t6)*((species0 * fbfp)-yb_t6), (float)0.5);
		yb_s += (float)pow((float)(tstep == t8)*((species0 * fbfp)-yb_t8)*((species0 * fbfp)-yb_t8), (float)0.5);
		yb_s += (float)pow((float)(tstep == t10)*((species0 * fbfp)-yb_t10)*((species0 * fbfp)-yb_t10), (float)0.5);
		

		ym_s += (float)pow((float)(tstep == t2)*((species1 * fmct)-ym_t2)*((species1 * fmct)-ym_t2), (float)0.5);
		ym_s += (float)pow((float)(tstep == t4)*((species1 * fmct)-ym_t4)*((species1 * fmct)-ym_t4), (float)0.5);
		ym_s += (float)pow((float)(tstep == t6)*((species1 * fmct)-ym_t6)*((species1 * fmct)-ym_t6), (float)0.5);
		ym_s += (float)pow((float)(tstep == t8)*((species1 * fmct)-ym_t8)*((species1 * fmct)-ym_t8), (float)0.5);
		ym_s += (float)pow((float)(tstep == t10)*((species1 * fmct)-ym_t10)*((species1 * fmct)-ym_t10), (float)0.5);
		


	}
	
	ELEMENT[BLINDEX].errorbfp[THINDEX] = yb_s;
	ELEMENT[BLINDEX].errormct[THINDEX] = ym_s;
	
	
}


/*
	Combined Taustep random number generator
	GPU Gems 3: http://http.developer.nvidia.com/GPUGems3/gpugems3_ch37.html
*/
__device__ unsigned TausStep(unsigned *z, int S1, int S2, int S3, unsigned M)  
{  
  unsigned b=(((*z << S1) ^ *z) >> S2);  
  return *z =  (((*z & M) << S3) ^ b);
}  

__device__ unsigned LCGStep(unsigned *z, unsigned A, unsigned C)  
{  
	return (A*(*z)+C);  
}  

__device__ float HybridTaus(unsigned *z1, unsigned *z2, unsigned *z3, unsigned *z4)  
{  
  return 2.3283064365387e-10 * (              // Periods  
    TausStep(z1, 13, 19, 12, 4294967294UL) ^  // p1=2^31-1  
    TausStep(z2, 2, 25, 4, 4294967288UL) ^    // p2=2^30-1  
    TausStep(z3, 3, 11, 17, 4294967280UL) ^   // p3=2^28-1  
    LCGStep(z4, 1664525, 1013904223UL)        // p4=2^32  
   );  
} 




/**************************************
 * 	Triple induction one experiment fit
 **************************************/
__global__ void ssaKernel9(curandState * STATE, heapElement * ELEMENT, int seed){

	//Register
	#include "exp8_parallel_imp\ssaspecies.h"

	//Random and indexing;
	int THINDEX = threadIdx.x;
	int BLINDEX = blockIdx.x;
	
	//Parameters
	float ta				= 1.0*ELEMENT[BLINDEX].vector[0];					//binding probability
	float tda				= 0.0001*ELEMENT[BLINDEX].vector[1];				//unbinding probability
	float tb				= 1.0*ELEMENT[BLINDEX].vector[2];					//binding probability
	float tdb				= 0.0001*ELEMENT[BLINDEX].vector[3];				//unbinding probability
	float tdg				= 100*0.00001*ELEMENT[BLINDEX].vector[4];			//degradation probability
	float r_vp16_vp16		= 0.0001*ELEMENT[BLINDEX].vector[5];				//fully activated transcription
	float r_vp16_krab		= 0.01*r_vp16_vp16*ELEMENT[BLINDEX].vector[12];		//competitive binding
	float tk				= 4.0*ELEMENT[BLINDEX].vector[14];					//KRAB right of way
	float rdg				= 0.00001*ELEMENT[BLINDEX].vector[6];				//reporter degradation
	float cta				= 0.5*r_vp16_vp16*0.1*ELEMENT[BLINDEX].vector[10];	//constitutive promoter transcription, LOWER FREQUENCY
	float p_leak_pir		= 10.0*0.1*ELEMENT[BLINDEX].vector[7];				//pir promotor leakage
	float p_leak_etr		= 10.0*0.1*ELEMENT[BLINDEX].vector[8];				//etr promotor leakage
	float tas				= 1.0*ELEMENT[BLINDEX].vector[11];					//TAL A/B asymmetry
	float pas				= 1.0*ELEMENT[BLINDEX].vector[13];					//Induction asymmetry
	float indmax			= 100.0*ELEMENT[BLINDEX].vector[16];				//Induction

	//Float
	float fbfp = 0.001*ELEMENT[BLINDEX].vector[9];
	float fmct = 0.001*ELEMENT[BLINDEX].vector[9];
	
	//Sample points
	//WARNING! Hard coding.
	int day		  = 18000;
	int t1		= day*3; int t2 = day*6; int t3 = day*9; int t4 = day*13; int t5 = day*17; 
	float yb_s	= 0;
	float ym_s	= 0;
	
	//Induction
	//int induced = 1 + 1*(THINDEX > HISTSIZE/2.0);
	int induced = 1;

	//Target values

	//TALB:KRAB - BFP
	float yb_t1 = 0.274;
	float yb_t2 = 0.086;
	float yb_t3 = 0.051;
	float yb_t4 = 0.149;
	float yb_t5 = 0.094;


	//TALA:KRAB - MCT
	float ym_t1 = 0.025;
	float ym_t2 = 0.013;
	float ym_t3 = 0.024;
	float ym_t4 = 0.115;
	float ym_t5 = 0.067;

	
	//Load vectors
	#include "exp8_parallel_imp\ssacoefs.h"
	
	
	//High log
	#define HL_AGM 20
	#define HL_PI 3.14159
	#define HL_LN2 0.69314
	#define HL_MBITS 8
	float hl_f;
	float hl_a;
	float hl_g;
	float hl_at;
	float hl_s;
	float hl_x;
	float hl_y;
	float hl_lnf;

		
	//Misc
	//int i		= 0;
	int j		= 0;
	int mi		= 0;
	int tstep	= 0;
	float tau	= 0;
	float t	= 0;
	
	//Ssa misc
	float r1;
	float r2;
	float wsum = 0;
	#include "exp8_parallel_imp\ssawthreshinit.h"
	float wtcumsum = 0;

	//Random number generator;
	unsigned z1;
	unsigned z2;
	unsigned z3;
	unsigned z4;

	z1 = 128+seed*THINDEX*BLINDEX;
	z2 = 129+seed*THINDEX*BLINDEX;
	z3 = 130+seed*THINDEX*BLINDEX;
	z4 = seed*THINDEX*BLINDEX;

	while(tstep <= STEPCOUNT){

		r1 = HybridTaus(&z1,&z2,&z3,&z4);
		r2 = HybridTaus(&z1,&z2,&z3,&z4);

		//Include propensities calculations
		#include "exp8_parallel_imp\ssaweights.h"
	
		//Calculate wsum
		//weight thresholds
		wsum = 0;
		wtcumsum = 0;
		wsum = wsum  + w0 + w1 + w2 + w3 + w4 + w5 + w6 + w7 + w8 + w9 + w10 + w11 + w12 + w13 + w14 + w15 + w16 + w17 + w18 + w19 + w20 + w21 + w22 + w23 + w24 + w25 + w26 + w27 + w28 + w29 + w30 + w31 + w32 + w33 + w34 + w35 + w36 + w37 + w38 + w39 + w40 + w41;
		#include "exp8_parallel_imp\ssawthresh.h"


		////Determine tau;
		//tau = 1.0/wsum * logf(1.0/(1-r1));
		hl_f = 1/(1-r1);
		hl_s = hl_f*pow(2.0,HL_MBITS);		
		hl_y = 4.0/hl_s;
		hl_x = 1;
		//agm start
		hl_a = 0.5*(hl_x+hl_y);
		hl_g = sqrt(hl_x*hl_y);
		for(j=0; j<HL_AGM; j++){
			hl_at = hl_a;
			hl_a = 0.5*(hl_a+hl_g);
			hl_g = sqrt(hl_g*hl_at);
		}
		//agm end
		hl_lnf = HL_PI / (2*hl_a) - HL_MBITS * HL_LN2;
		tau = 1/wsum * hl_lnf;

		////Determine reaction index.
		mi = -1;
		#include "exp8_parallel_imp\ssawtmi.h"
		

		////Update species matrix;
		#include "exp8_parallel_imp\ssadifmx.h"
		

		//Update time and step
		t = t+tau;
		tstep++;

		//switching
		//1: er_pi
		//2: pi_er
		species26 = (induced == 1 && tstep < t2)*indmax;
		species27 = (induced == 2 && tstep < t2)*indmax;

		//Sample values
		yb_s += (float)pow((float)(tstep == t1)*((species0 * fbfp)-yb_t1)*((species0 * fbfp)-yb_t1), (float)1.0); 
		yb_s += (float)pow((float)(tstep == t2)*((species0 * fbfp)-yb_t2)*((species0 * fbfp)-yb_t2), (float)1.0);
		yb_s += (float)pow((float)(tstep == t3)*((species0 * fbfp)-yb_t3)*((species0 * fbfp)-yb_t3), (float)1.0);
		yb_s += (float)pow((float)(tstep == t4)*((species0 * fbfp)-yb_t4)*((species0 * fbfp)-yb_t4), (float)1.0);
		yb_s += (float)pow((float)(tstep == t5)*((species0 * fbfp)-yb_t5)*((species0 * fbfp)-yb_t5), (float)1.0);
		

		ym_s += (float)pow((float)(tstep == t1)*((species1 * fmct)-ym_t1)*((species1 * fmct)-ym_t1), (float)1.0);
		ym_s += (float)pow((float)(tstep == t2)*((species1 * fmct)-ym_t2)*((species1 * fmct)-ym_t2), (float)1.0);
		ym_s += (float)pow((float)(tstep == t3)*((species1 * fmct)-ym_t3)*((species1 * fmct)-ym_t3), (float)1.0);
		ym_s += (float)pow((float)(tstep == t4)*((species1 * fmct)-ym_t4)*((species1 * fmct)-ym_t4), (float)1.0);
		ym_s += (float)pow((float)(tstep == t5)*((species1 * fmct)-ym_t5)*((species1 * fmct)-ym_t5), (float)1.0);
		


	}
	
	ELEMENT[BLINDEX].errorbfp[THINDEX] = yb_s;
	ELEMENT[BLINDEX].errormct[THINDEX] = ym_s;
	
	
}